This document contains a description of optimisation tools for electric system planing simulation. If you are only interested in operation you might prefer to go back to file optim-Operation.ipynb. In files BasicFunctionalities/input-XXX.ipynb you can learn to understand input data (consumption, availability).
(Text below is almost the same as for operation)
This document will gives a chance to understand
It proposes to enter the subject by increasing progressively the number of variables and constraints in the optimisation problem, hence moving toward more realism through the document, introducing:
It relies on different test cases that allow to
If, after reading this file, you want to build your own pyomo model you can go to optim-Planing-Advanced.ipynb.
#region importation of modules
import os
if os.path.basename(os.getcwd())=="BasicFunctionalities":
os.chdir('..') ## to work at project root like in any IDE
import sys
if sys.platform != 'win32':
myhost = os.uname()[1]
else : myhost = ""
if (myhost=="jupyter-sop"):
## for https://jupyter-sop.mines-paristech.fr/ users, you need to
# (1) run the following in a terminal
if (os.system("/opt/mosek/9.2/tools/platform/linux64x86/bin/lmgrd -c /opt/mosek/9.2/tools/platform/linux64x86/bin/mosek.lic -l lmgrd.log")==0):
os.system("/opt/mosek/9.2/tools/platform/linux64x86/bin/lmutil lmstat -c 27007@127.0.0.1 -a")
# (2) definition of license
os.environ["MOSEKLM_LICENSE_FILE"] = '@jupyter-sop'
import numpy as np
import pandas as pd
import csv
#import docplex
import datetime
import copy
import plotly.graph_objects as go
import matplotlib.pyplot as plt
from sklearn import linear_model
import sys
from functions.f_planingModels import *
from functions.f_optimization import *
from functions.f_graphicalTools import *
from functions.f_consumptionModels import *
# Change this if you have other solvers obtained here
## https://ampl.com/products/solvers/open-source/
## for eduction this site provides also several professional solvers, that are more efficient than e.g. cbc
#endregion
#region Solver and data location definition
InputFolder='Data/input/'
solver= 'mosek' ## no need for solverpath with mosek.
BaseSolverPath='/Users/robin.girard/Documents/Code/Packages/solvers/ampl_macosx64'
sys.path.append(BaseSolverPath)
solvers= ['gurobi','knitro','cbc'] # 'glpk' is too slow 'cplex' and 'xpress' do not work
solverpath= {}
for solver in solvers : solverpath[solver]=BaseSolverPath+'/'+solver
cplexPATH='/Applications/CPLEX_Studio1210/cplex/bin/x86-64_osx'
sys.path.append(cplexPATH)
solverpath['cplex']=cplexPATH+"/"+"cplex"
solver = 'mosek'
#endregion
Before you start with the math, you should
#region I - Simple single area : loading parameters
Zones="FR" ; year=2013
#### reading areaConsumption availabilityFactor and TechParameters CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Planing-Simple_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
#### Selection of subset
Selected_TECHNOLOGIES=['OldNuke','CCG'] #you can add technologies here
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
#TechParameters.loc[TechParameters.TECHNOLOGIES=="OldNuke",'maxCapacity']=63000 ## Limit to actual installed capacity
#endregion
#region I - Simple single area : Solving and loading results
model = GetElectricSystemModel_PlaningSingleNode(areaConsumption,availabilityFactor,TechParameters)
if solver in solverpath : opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
## result analysis
Variables=getVariables_panda_indexed(model)
print(extractCosts(Variables))
print(extractEnergyCapacity(Variables))
#pour avoir la production en KWh de chaque moyen de prod chaque heure
production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
### Check sum Prod = Consumption
Delta=(production_df.sum(axis=1) - areaConsumption.areaConsumption);
abs(Delta).max()
#endregion
Capacity_Milliards_euros Energy_Milliards_euros
TECHNOLOGIES
OldNuke 19.060436 3.308393
CCG 2.628352 1.333826
Capacity_GW Production_TWh
TECHNOLOGIES
OldNuke 77.778651 472.627621
CCG 22.485684 19.386999
0.0
Verify that the sum of market prices allows all actors to cover fixed and marginal cost. do they earn more ? why ?
#region I - Simple single area : visualisation and lagrange multipliers
### representation des résultats
fig=MyStackedPlotly(y_df=production_df,Conso = areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#### lagrange multipliers
Constraints= getConstraintsDual_panda(model)
# Analyse energyCtr
energyCtrDual=Constraints['energyCtr']; energyCtrDual['energyCtr']=energyCtrDual['energyCtr']
energyCtrDual
round(energyCtrDual.energyCtr,2).unique()
# Analyse CapacityCtr
CapacityCtrDual=Constraints['CapacityCtr'].pivot(index="Date",columns='TECHNOLOGIES', values='CapacityCtr')*1000000;
round(CapacityCtrDual,2)
round(CapacityCtrDual.OldNuke,2).unique() ## if you increase by Delta the installed capacity of nuke you decrease by xxx the cost when nuke is not sufficient
round(CapacityCtrDual.CCG,2).unique() ## increasing the capacity of CCG as no effect on prices
#endregion
array([-0.0000e+00, -1.1689e+11])
In the this section, we will increase the complexity of the problem given in Section 2 and add : dependency on area z (country),a congestion constraint, ramp constraints.
\begin{align} &\text{Cost function }& &\min_{x} \sum_z \left ( \beta_{iz}\bar{x_{iz}} + \sum_t \sum_i \pi_{iz} x_{itz} \right ) \;\;\; & & \pi_{iz} \text{ marginal cost}\\ &\text{Power limit } & &\text{ s.t.} \;\; 0 \leq x_{itz}\leq a_{itz} \bar{x_{iz}} & &\bar{x_{iz}} \text{ installed power, } a_{itz} \text{ availability}\\ &\text{Meet demand } & & \sum_i x_{itz} \geq C_{tz} && C_{tz} \text{ Consumption}\\ &\text{Stock limit } & &\sum_t x_{it}\leq E_i && E_i=\bar{x_i}*N_i \text{ Energy capacity limit}\\ &\text{ramp limit } & &rc^-_i *x_{it}\leq x_{it}-x_{i(t+1)}\leq rc^+_i *x_{it} && rc^+_i rc^-_i\text{ ramp limit}\\ \end{align}#region II - Ramp Single area : loading parameters loading parameterscase with ramp constraints
Zones="FR"
year=2013
Selected_TECHNOLOGIES=['OldNuke', 'CCG',"curtailment"] #you'll add 'Solar' after
#### reading CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Planing-RAMP1BIS_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
TechParameters.loc["OldNuke",'RampConstraintMoins']=0.01 ## a bit strong to put in light the effect
TechParameters.loc["OldNuke",'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
#endregion
#region II - Ramp Single area : solving and loading results
model = GetElectricSystemModel_PlaningSingleNode(areaConsumption,availabilityFactor,TechParameters)
opt = SolverFactory(solver)
results=opt.solve(model)
Variables=getVariables_panda_indexed(model)
print(extractCosts(Variables))
print(extractEnergyCapacity(Variables))
#pour avoir la production en KWh de chaque moyen de prod chaque heure
production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
### Check sum Prod = Consumption
Delta=(production_df.sum(axis=1) - areaConsumption.areaConsumption);
abs(Delta).max()
#endregion
Capacity_Milliards_euros Energy_Milliards_euros
TECHNOLOGIES
OldNuke 15.43878 2.701839
curtailment 0.00000 910.235379
CCG 1.40268 5.207915
Capacity_GW Production_TWh
TECHNOLOGIES
OldNuke 63.0 385.977000
curtailment 1000.0 30.341179
CCG 12.0 75.696441
1.4551915228366852e-11
#region II - Ramp Single area : visualisation and lagrange multipliers
fig=MyStackedPlotly(y_df=production_df,Conso = areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#### lagrange multipliers
Constraints= getConstraintsDual_panda(model)
# Analyse energyCtr
energyCtrDual=Constraints['energyCtr']; energyCtrDual['energyCtr']=energyCtrDual['energyCtr']*1000000
energyCtrDual
round(energyCtrDual.energyCtr,2).unique()
# Analyse CapacityCtr
CapacityCtrDual=Constraints['CapacityCtr'].pivot(index="Date",columns='TECHNOLOGIES', values='CapacityCtr')*1000000;
round(CapacityCtrDual,2)
round(CapacityCtrDual.OldNuke,2).unique() ## if you increase by Delta the installed capacity of nuke you decrease by xxx the cost when nuke is not sufficient
round(CapacityCtrDual.CCG,2).unique() ## increasing the capacity of CCG as no effect on prices
#endregion
array([-0.00000e+00, -2.63468e+10, -2.99312e+10, -2.98694e+10,
-2.87570e+10, -2.90660e+10, -2.89424e+10, -2.98076e+10,
-2.93132e+10, -2.86334e+10, -2.85098e+10, -2.83862e+10,
-2.96840e+10, -2.97458e+10, -2.91278e+10, -2.95604e+10,
-2.93750e+10, -2.92514e+10, -2.94368e+10, -2.69648e+10,
-2.96222e+10, -2.91896e+10, -2.94986e+10, -2.59142e+10,
-2.88806e+10, -2.69030e+10, -2.71502e+10, -2.36276e+10,
-2.57906e+10, -2.88188e+10, -2.86952e+10, -3.48080e+09,
-2.85716e+10, -2.42456e+10, -2.25152e+10, -2.56670e+10,
-2.84480e+10, -2.28860e+10, -2.28242e+10, -1.51610e+10,
-2.56052e+10, -2.30714e+10, -2.58524e+10, -2.90042e+10,
-2.82626e+10])
Math here are in 3.1
#region III - Ramp multiple area : loading parameters
Zones="FR_DE_GB_ES"
year=2016
Selected_AREAS=["FR","DE"]
Selected_TECHNOLOGIES=['OldNuke', 'CCG','WindOnShore',"curtailment"] #you'll add 'Solar' after #'NewNuke', 'HydroRiver', 'HydroReservoir','WindOnShore', 'WindOffShore', 'Solar', 'Curtailement'}
#### reading CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["AREAS","Date"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["AREAS","Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Planing_MultiNode_DE-FR_TECHNOLOGIES_AREAS.csv',
sep=',',decimal='.',skiprows=0,comment="#").set_index(["AREAS","TECHNOLOGIES"])
ExchangeParameters = pd.read_csv(InputFolder+'Hypothese_DE-FR_AREAS_AREAS.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["AREAS","AREAS.1"])
#### Selection of subset
TechParameters=TechParameters.loc[(Selected_AREAS,Selected_TECHNOLOGIES),:]
areaConsumption=areaConsumption.loc[(Selected_AREAS,slice(None)),:]
availabilityFactor=availabilityFactor.loc[(Selected_AREAS,slice(None),Selected_TECHNOLOGIES),:]
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintMoins']=0.01 ## a bit strong to put in light the effect
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
#endregion
#region III - Ramp multiple area : solving and loading results
### small data cleaning
availabilityFactor.availabilityFactor[availabilityFactor.availabilityFactor>1]=1
model = GetElectricSystemModel_PlaningMultiNode(areaConsumption,availabilityFactor,TechParameters,ExchangeParameters)
opt = SolverFactory(solver)
results=opt.solve(model)
Variables=getVariables_panda_indexed(model)
print(extractCosts(Variables))
#print(extractEnergyCapacity(Variables)) #bug avec cette ligne à corriger
production_df=EnergyAndExchange2Prod(Variables)
### Check sum Prod = Consumption
Delta= production_df.sum(axis=1)-areaConsumption.areaConsumption
abs(Delta).sum()
## adding dates
fig=MyAreaStackedPlot(production_df,Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
Constraints= getConstraintsDual_panda(model)
Constraints.keys()
#endregion
Capacity_Milliards_euros Energy_Milliards_euros
AREAS TECHNOLOGIES
DE OldNuke 15.476304 2.584650
curtailment 0.000000 0.074833
CCG 5.041756 3.556250
WindOnShore 6.360066 0.000000
FR OldNuke 17.845299 3.026397
curtailment 0.000000 1.118105
CCG 2.503603 1.993053
WindOnShore 1.784651 0.000000
dict_keys(['energyCostsDef', 'capacityCostsCtr', 'CapacityCtr', 'exchangeCtrPlus', 'exchangeCtrMoins', 'exchangeCtr2', 'energyCtr', 'maxCapacityCtr', 'minCapacityCtr', 'storageCtr', 'rampCtrPlus', 'rampCtrMoins', 'rampCtrPlus2', 'rampCtrMoins2'])
Just have a look at optim-Storage.ipynb.
You can find the maths in optim-Operation.ipynb. This only change is that this time, storage parameters $P_{max}$ and $C_{max}$ are decision variables.
#region IV Ramp+Storage single area : loading parameters
Zones="FR"
year=2013
Selected_TECHNOLOGIES=['OldNuke','WindOnShore', 'CCG',"curtailment"] ## try adding 'HydroRiver', 'HydroReservoir'
#### reading CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Planing-RAMP1BIS_TECHNOLOGIES.csv',
sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
StorageParameters = pd.read_csv(InputFolder+'Planing-RAMP1_STOCK_TECHNO.csv',sep=',',decimal='.',skiprows=0).set_index(["STOCK_TECHNO"])
#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
#TechParameters.loc["CCG",'capacity']=100000 ## margin to make everything work
TechParameters.loc["OldNuke",'RampConstraintMoins']=0.02 ## a bit strong to put in light the effect
TechParameters.loc["OldNuke",'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
#endregion
#region IV Ramp+Storage single area : solving and loading results
model= GetElectricSystemModel_PlaningSingleNode_withStorage(areaConsumption,availabilityFactor,
TechParameters,StorageParameters)
if solver in solverpath : opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
Variables = getVariables_panda_indexed(model)
Constraints = getConstraintsDual_panda(model)
production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
production_df.loc[:,'Storage'] = Variables['storageOut'].pivot(index='Date',columns='STOCK_TECHNO',values='storageOut').sum(axis=1)-Variables['storageIn'].pivot(index='Date',columns='STOCK_TECHNO',values='storageIn').sum(axis=1) ### put storage in the production time series
production_df.sum(axis=0)/10**6 ### energies produites TWh
production_df[production_df>0].sum(axis=0)/10**6 ### energies produites TWh
production_df.max(axis=0)/1000 ### Pmax en GW
fig=MyStackedPlotly(y_df=production_df, Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#endregion
Case V is in section 5. (lots of renewable)
#region VI Case Storage + CCG + Nuke (Ramp+Storage single area) : loading parameters
Zones="FR"
year=2013
Selected_TECHNOLOGIES=['CCG', 'NewNuke',"curtailment",'HydroRiver', 'HydroReservoir']
#### reading CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Planing-RAMP1BIS_TECHNOLOGIES.csv',
sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
StorageParameters = pd.read_csv(InputFolder+'Planing-RAMP1_STOCK_TECHNO.csv',sep=',',decimal='.',skiprows=0).set_index(["STOCK_TECHNO"])
#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
TechParameters.loc["CCG",'energyCost']=300
TechParameters.loc["CCG",'RampConstraintMoins']=0.05 ## a bit strong to put in light the effect
TechParameters.loc["CCG",'RampConstraintPlus']=0.05 ## a bit strong to put in light the effect
TechParameters.loc["NewNuke",'RampConstraintMoins']=0.01 ## a bit strong to put in light the effect
TechParameters.loc["NewNuke",'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
#endregion
#region VI Case Storage + CCG + Nuke (Ramp+Storage single area) : solving and loading results
model= GetElectricSystemModel_PlaningSingleNode_withStorage(areaConsumption,availabilityFactor,
TechParameters,StorageParameters)
if solver in solverpath : opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
Variables = getVariables_panda_indexed(model)
Constraints = getConstraintsDual_panda(model)
print(extractCosts(Variables))
print(extractEnergyCapacity(Variables))
production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
stockage=Variables['storageOut'].pivot(index='Date',columns='STOCK_TECHNO',values='storageOut').sum(axis=1)-Variables['storageIn'].pivot(index='Date',columns='STOCK_TECHNO',values='storageIn').sum(axis=1)
areaConsumption['Storage']=stockage
areaConsumption['NewConsumption']=areaConsumption['areaConsumption']-stockage
Delta= production_df.sum(axis=1)-areaConsumption["NewConsumption"]
print(sum(abs(Delta)))
production_df.loc[:,'Storage'] = areaConsumption["Storage"] ### put storage in the production time series
production_df.sum(axis=0)/10**6 ### energies produites TWh
production_df[production_df>0].sum(axis=0)/10**6 ### energies produites TWh
production_df.max(axis=0)/1000 ### Pmax en GW
fig=MyStackedPlotly(y_df=production_df, Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#endregion
Capacity_Milliards_euros Energy_Milliards_euros
TECHNOLOGIES
HydroReservoir 0.565250 0.000000
curtailment 0.000000 0.002137
CCG 1.122162 1.087435
NewNuke 27.966525 2.986631
HydroRiver 0.807500 0.000000
Capacity_GW Production_TWh
TECHNOLOGIES
HydroReservoir 7.000000 14.700000
curtailment 1000.000000 0.000071
CCG 9.600153 3.624784
NewNuke 70.903645 426.661526
HydroRiver 10.000000 48.798573
2.2922904463484883e-08
#region VI Ramp+Storage Multi area : loading parameters
Zones="FR_DE_GB_ES"
year=2016
Selected_AREAS=["FR","DE"]
Selected_TECHNOLOGIES=['OldNuke', 'CCG','WindOnShore',"curtailment"] #you'll add 'Solar' after #'NewNuke', 'HydroRiver', 'HydroReservoir','WindOnShore', 'WindOffShore', 'Solar', 'Curtailement'}
#### reading CSV files
TechParameters = pd.read_csv(InputFolder+'Planing_MultiNode_DE-FR_TECHNOLOGIES_AREAS.csv',
sep=',',decimal='.',comment="#").set_index(["AREAS","TECHNOLOGIES"])
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["AREAS","Date"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["AREAS","Date","TECHNOLOGIES"])
ExchangeParameters = pd.read_csv(InputFolder+'Hypothese_DE-FR_AREAS_AREAS.csv',
sep=',',decimal='.',skiprows=0,comment="#").set_index(["AREAS","AREAS.1"])
StorageParameters = pd.read_csv(InputFolder+'Planing_MultiNode_AREAS_DE-FR_STOCK_TECHNO.csv',sep=',',decimal='.',comment="#",skiprows=0).set_index(["AREAS","STOCK_TECHNO"])
#### Selection of subset
TechParameters=TechParameters.loc[(Selected_AREAS,Selected_TECHNOLOGIES),:]
areaConsumption=areaConsumption.loc[(Selected_AREAS,slice(None)),:]
availabilityFactor=availabilityFactor.loc[(Selected_AREAS,slice(None),Selected_TECHNOLOGIES),:]
TechParameters.loc[(slice(None),'CCG'),'energyCost']=300 ## margin to make everything work
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintMoins']=0.01 ## a bit strong to put in light the effect
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
#endregion
#region VI Ramp+Storage multi area : solving and loading results
model= GetElectricSystemModel_PlaningMultiNode_withStorage(areaConsumption,availabilityFactor,
TechParameters,ExchangeParameters,StorageParameters)
if solver in solverpath : opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
Variables = getVariables_panda_indexed(model)
Constraints = getConstraintsDual_panda(model)
production_df=EnergyAndExchange2Prod(Variables)
stockage=Variables['storageOut'].pivot(index=['AREAS','Date'],columns='STOCK_TECHNO',values='storageOut').sum(axis=1)-Variables['storageIn'].pivot(index=['AREAS','Date'],columns='STOCK_TECHNO',values='storageIn').sum(axis=1)
areaConsumption['Storage']=stockage
areaConsumption['NewConsumption']=areaConsumption['areaConsumption']-stockage
Delta=(production_df.sum(axis=1) - areaConsumption.NewConsumption); ## comparaison à la conso incluant le stockage
abs(Delta).max()
production_df.loc[:,'Storage'] = -areaConsumption["Storage"] #### ajout du stockage comme production
Delta=(production_df.sum(axis=1) - areaConsumption.areaConsumption);
abs(Delta).max()
fig=MyAreaStackedPlot(production_df,Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
abs(areaConsumption["Storage"]).groupby(by="AREAS").sum() ## stockage
production_df.groupby(by="AREAS").sum()/10**6 ### energies produites TWh
production_df[production_df>0].groupby(by="AREAS").sum()/10**6 ### energies produites TWh
production_df.groupby(by="AREAS").max()/1000 ### Pmax en GW ### le stockage ne fait rien en Allemagne ??? bizarre
production_df.groupby(by="AREAS").min()/1000 ### Pmax en GW
#endregion
| CCG | OldNuke | WindOnShore | curtailment | DE | FR | Storage | |
|---|---|---|---|---|---|---|---|
| AREAS | |||||||
| DE | -1.017469e-16 | 14.032457 | 0.0 | -0.0 | 0.0 | -9.0 | -5.0 |
| FR | -4.747401e-17 | 27.833874 | 0.0 | -0.0 | -9.0 | 0.0 | -5.0 |
#region V Case Storage + CCG + PV + Wind + hydro (Ramp+Storage single area) : loading parameters
Zones="FR"
year=2013
Selected_TECHNOLOGIES=['CCG', 'WindOnShore','WindOffShore','Solar',"curtailment",'HydroRiver', 'HydroReservoir']
#### reading CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Planing-RAMP1BIS_TECHNOLOGIES.csv',
sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
StorageParameters = pd.read_csv(InputFolder+'Planing-RAMP1_STOCK_TECHNO.csv',sep=',',decimal='.',skiprows=0).set_index(["STOCK_TECHNO"])
#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
TechParameters.loc["CCG",'maxCapacity']=100000 ## margin to make everything work
TechParameters.loc["CCG",'energyCost']=300
TechParameters.loc["CCG",'RampConstraintMoins']=0.05 ## a bit strong to put in light the effect
TechParameters.loc["CCG",'RampConstraintPlus']=0.05 ## a bit strong to put in light the effect
#endregion
#region V Case Storage + CCG + PV + Wind (Ramp+Storage single area) : solving and loading results
model= GetElectricSystemModel_PlaningSingleNode_withStorage(areaConsumption,availabilityFactor,
TechParameters,StorageParameters)
if solver in solverpath : opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
Variables = getVariables_panda_indexed(model)
Constraints = getConstraintsDual_panda(model)
production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
production_df.loc[:,'Storage'] = Variables['storageOut'].pivot(index='Date',columns='STOCK_TECHNO',values='storageOut').sum(axis=1)-Variables['storageIn'].pivot(index='Date',columns='STOCK_TECHNO',values='storageIn').sum(axis=1) ### put storage in the production time series
#production_df.sum(axis=0)/10**6 ### energies produites TWh
print(production_df[production_df>0].sum(axis=0)/10**6) ### energies produites TWh
print(production_df.max(axis=0)/1000) ### Pmax en GW
fig=MyStackedPlotly(y_df=production_df, Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#endregion
TECHNOLOGIES CCG 57.195768 HydroReservoir 14.700000 HydroRiver 45.555465 Solar 35.627231 WindOffShore 0.000000 WindOnShore 342.266979 curtailment 0.004264 Storage 7.854631 dtype: float64 TECHNOLOGIES CCG 57.749251 HydroReservoir 7.000000 HydroRiver 8.291000 Solar 47.076554 WindOffShore 0.000000 WindOnShore 84.071000 curtailment 2.490076 Storage 5.000000 dtype: float64
#region VII - Simple single area +4 million EV + demande side management +30TWh H2: loading parameters
Zones="FR" ; year=2013
#### reading areaConsumption availabilityFactor and TechParameters CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["Date","TECHNOLOGIES"])
#obtaining industry-metal consumption
Profile_df=pd.read_csv(InputFolder+"ConsumptionDetailedProfiles.csv").set_index(["Mois", "heures",'Nature', 'type', 'UsagesGroupe', 'UsageDetail', "WeekDay"])
Profile_df_merged=ComplexProfile2Consumption(Profile_df,areaConsumption,TemperatureName='areaConsumption') #TODO : change ComplexProfile2Consumption to adapt to the new date format
#Profile_df_merged_spread_0=Profile_df_merged.groupby(["Date","type"]).sum().reset_index().drop(columns=["areaConsumption"]).pivot(index="Date", columns=['type'], values='Conso');
Profile_df_merged_spread = Profile_df_merged.groupby(["Date","Nature","UsagesGroupe","type"]).sum().reset_index().drop(columns=["areaConsumption"]).pivot(index="Date", columns=["Nature",'type',"UsagesGroupe"], values='Conso');
steel_consumption=Profile_df_merged_spread.loc[:,("MineraiMetal","Ind",'Process')]*areaConsumption.loc[:,"areaConsumption"]
steel_consumption.iloc[0]=110
steel_consumption.max()
steel_consumption[steel_consumption.isna()]=110
steel_consumption.isna().sum()
# if you want to change thermal sensitivity + add electric vehicle
ConsoTempe_df=pd.read_csv(InputFolder+'ConsumptionTemperature_1996TO2019_FR.csv',parse_dates=['Date']).set_index(["Date"]) #
ConsoTempe_df_nodup=ConsoTempe_df.loc[~ConsoTempe_df.index.duplicated(),:]
(ConsoTempeYear_decomposed_df,Thermosensibilite)=Decomposeconso(areaConsumption.join(ConsoTempe_df_nodup)[['areaConsumption','Temperature']],
TemperatureThreshold=15,ConsumptionName="areaConsumption")
areaConsumption = areaConsumption.assign(areaConsumption = ConsoTempeYear_decomposed_df.loc[:,'TS_C']+ConsoTempeYear_decomposed_df.loc[:,'NTS_C'])#+VE_consumption.loc[:,'Consumption'])
VEProfile_df=pd.read_csv(InputFolder+'EVModel.csv', sep=';')
NbVE=10 # millions
ev_consumption = NbVE*Profile2Consumption(Profile_df=VEProfile_df,Temperature_df = ConsoTempe_df_nodup.loc[str(year)][['Temperature']])[['Consumption']]
TechParameters = pd.read_csv(InputFolder+'Planing-RAMP1BIS_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["TECHNOLOGIES"])
StorageParameters = pd.read_csv(InputFolder + 'Planing-RAMP1_STOCK_TECHNO.csv', sep=',', decimal='.',
skiprows=0).set_index(["STOCK_TECHNO"])
ConsoParameters = pd.read_csv(InputFolder + "Planing-Conso-FLEX_CONSUM.csv", sep=";").set_index(["FLEX_CONSUM"])
Selected_TECHNOLOGIES=['OldNuke','CCG','TAC', 'WindOnShore', 'WindOffShore','HydroReservoir','HydroRiver','Solar','curtailment']#you can add technologies here
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
h2_Energy = 30000## H2 volume in GWh/year
h2_Energy_flat_consumption = ev_consumption.Consumption*0+h2_Energy/8760
to_flexible_consumption=pd.DataFrame({'areaConsumption': steel_consumption,'FLEX_CONSUM' : 'Steel'}).reset_index().set_index(['Date','FLEX_CONSUM']).\
append(pd.DataFrame({'areaConsumption': ev_consumption.Consumption,'FLEX_CONSUM' : 'EV'}).reset_index().set_index(['Date','FLEX_CONSUM'])).\
append(pd.DataFrame({'areaConsumption': h2_Energy_flat_consumption,'FLEX_CONSUM' : 'H2'}).reset_index().set_index(['Date','FLEX_CONSUM']))
# endregion
ConsoParameters.loc['H2','LoadCost']
# €/kW/an coût fixe additionnel pour un GW d'électrolyseur en plus en supposant que l'on construit
model= Model_SingleNode_online_flex(areaConsumption, availabilityFactor, to_flexible_consumption,
ConsoParameters, TechParameters, StorageParameters)
if solver in solverpath : opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
Variables = getVariables_panda_indexed(model)
Variables.keys()
Variables['increased_max_power'] ## on a ajouté X GW à ce qui existait.
print(extractCosts(Variables))
print(extractEnergyCapacity(Variables))
Constraints = getConstraintsDual_panda(model)
production_df=Variables['energy'].pivot(index="Date",columns='TECHNOLOGIES', values='energy')
production_df.loc[:,'Storage'] = Variables['storageOut'].pivot(index='Date',columns='STOCK_TECHNO',values='storageOut').sum(axis=1)-Variables['storageIn'].pivot(index='Date',columns='STOCK_TECHNO',values='storageIn').sum(axis=1) ### put storage in the production time series
#production_df.sum(axis=0)/10**6 ### energies produites TWh
print(production_df[production_df>0].sum(axis=0)/10**6) ### energies produites TWh
print(production_df.max(axis=0)/1000) ### Pmax en GW
fig=MyStackedPlotly(y_df=production_df, Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
Capacity_Milliards_euros Energy_Milliards_euros
TECHNOLOGIES
HydroReservoir 0.565250 0.000000
curtailment 0.000000 0.078899
CCG 0.764257 1.118879
WindOnShore 1.645675 0.000000
WindOffShore 0.000000 0.000000
HydroRiver 0.807500 0.000000
OldNuke 15.438780 2.770718
TAC 1.007315 0.336349
Solar 0.608947 0.000000
Capacity_GW Production_TWh
TECHNOLOGIES
HydroReservoir 7.000000 14.700000
curtailment 1000.000000 0.002630
CCG 6.538254 16.262783
WindOnShore 13.410000 27.381252
WindOffShore 0.000000 0.000000
HydroRiver 10.000000 48.799242
OldNuke 63.000000 395.816809
TAC 10.053040 4.004153
Solar 8.640000 10.195112
TECHNOLOGIES
CCG 16.262783
HydroReservoir 14.700000
HydroRiver 48.799242
OldNuke 395.816809
Solar 10.195112
TAC 4.004153
WindOffShore 0.000000
WindOnShore 27.381252
curtailment 0.002630
Storage 5.306440
dtype: float64
TECHNOLOGIES
CCG 6.538254
HydroReservoir 7.000000
HydroRiver 8.291000
OldNuke 59.552112
Solar 6.524562
TAC 10.053040
WindOffShore 0.000000
WindOnShore 10.821131
curtailment 1.295601
Storage 5.000000
dtype: float64